---
title: "Rising temperatures erode human sleep globally: sleep impact projection plots"
author: "Kelton Minor*, Andreas Bjerre-Nielsen, Sigga Svala Jonasdottir, Sune Lehmann, Nick Obradovich"
date: "2/4/2022"
output: html_document
---


##Downscaled sleep impact projections for NASA NEX bias corrected and statistically downscaled CMIP5 climate model output run under RCP8.5 and RCP4.5

```{r, include=FALSE}
#Load packages
library(ggthemes) #themes for plotting
library(sandwich) #robust se estimator
library(lmtest) #lm functions
library(ordinal) #ologit funcitons
library(memisc) #import stata data
library(ggmap) #make nice plots with ggplot mapping
library(ggplot2) #nice plots
library(Hmisc) #multipurpose package
library(RColorBrewer) #make nice colors
library(grid) #plotting package
library(MASS) #multipurpose package
library(gplots) #nice plotting
library(lfe) #run linear fixed effects models
library(stargazer) #output nice tables
library(lubridate) #deal with dates
library(zoo) #deal with dates
library(reshape2) #powerful reshape package
library(sp) #general spatial class package
library(gstat) #geostat package
library(spacetime) #store data in proper spatial/temporal class
library(raster) #convert to raster and vice versa
library(maptools) #plot and modify shapefiles
library(rgeos) #for use with maptools for shapefiles
library(rgdal) #for reading in shapefiles
library(RSAGA) #needed for spatial downscaling
library(RCurl) #also needed for spatial downscaling
library(dplyr) #for data frame manipulation functions
library(plyr) #for additional manipulation functions
library(data.table) #for data.table functions
library(foreign) #for reading in .dta files
library(parallel) #for parallel processing functions
library(foreach) #for plyr parallel
library(doMC) #run multicore processes
library(caTools) #functions for fast running mean and sd
library(ncdf4) #tools for netCDF packages
library(stringr) #string tools
#library(kfigr) #figure referencing for markdown
library(knitr) #knitting
library(pander) #pandering
library(tidyr)
library(rasterVis) # plotting rasters
library(ggExtra) # for marginal histograms
library(gridExtra) # for multiplots
library(gtable) # for plots
library(bit64) # long integers
library(chron)
library(lattice)
library(ggridges)
library(ggplot2)
library(sf)
```

```{r}
#set path to your parent directory
pathdir <- "~/Replication_Code/"
```

#Excess Sleep Loss RCP4.5
```{r}
#1.Compute and plot NASA SEDAC GPWv4 population-weighted avg. annual temperature-attributed EXCESS SLEEP LOSS due to suboptimal temperature for every model-by-year (2010,2050,2099) using NASA NEX BCSD CMIP5 (RCP4.5) spline sleep projections

#load sleep loss projection data for RCP4.5 (spline model prediction run on climate model time-series)
geocoords_labels<-fread(paste0(pathdir,"model_data/global_downscaled_sleeploss_projection_RCP4.5_popweights.csv"))
geocoords_labels[,V1:=NULL] #drop index col

#compute population weights using NASA SEDAC 2015 GPWv4 population count data
geocoords_labels[,pop_tot:=sum(pop_count)] 
geocoords_labels[,pop_wts:=pop_count/pop_tot]

#for each model, multiply sleep loss by pop weight fraction (note: need to specify climate model columns)
cols<-names(geocoords_labels)[3:65] #select model cols
#geocoords_labels
geocoords_labels[,(cols):=lapply(.SD, "*",pop_wts),.SDcols=cols]
#geocoords_labels

#sum weighted sleep loss values
avg_model_year <- geocoords_labels[, lapply(.SD, sum), .SDcols = cols]
#avg_model_year

#compute avg. for 2050 and 2099 across all vars
avg_model_year_long <- melt(avg_model_year)

#scale effect to sleep loss (in hours) per person
avg_model_year_long$transformed <- avg_model_year_long[,.(transformed = (value/60))]

#create year variable 
avg_model_year_long$year <- str_sub(avg_model_year_long$variable, -7, -4)

#create model variable 
avg_model_year_long$model <- str_sub(avg_model_year_long$variable, -17, -8)

#add 2010 baseline
avg_model_2010<-unique(avg_model_year_long[year=="2010",.(model,loss2010=transformed)])
avg_model_year_long[,.N]
avg_model_year_long<-merge(avg_model_year_long,avg_model_2010,all.x=T,by=c("model"))
avg_model_year_long[,.N]

setnames(avg_model_year_long, "transformed", "annual_sleep_loss_per_person") #rename y var

#compute total annual sleep loss by model-year
avg_model_year_long[year!=2010,annual_sleep_loss_per_person:=(loss2010 + annual_sleep_loss_per_person)]
#avg_model_year_long
```

```{r}
#inspect min, median and max sleep loss projection values by year
avg_model_year_long[year=="2010",.(med=median(annual_sleep_loss_per_person), low = min(annual_sleep_loss_per_person), high = max(annual_sleep_loss_per_person))]
avg_model_year_long[year=="2050",.(med=median(annual_sleep_loss_per_person), low = min(annual_sleep_loss_per_person), high = max(annual_sleep_loss_per_person))]
avg_model_year_long[year=="2099",.(med=median(annual_sleep_loss_per_person), low = min(annual_sleep_loss_per_person), high = max(annual_sleep_loss_per_person))]
```

```{r}
#RCP4.5
#Plot population-weighted average sleep loss due to temperature (hours of sleep lost per person per year)
#average loss by year across all models (ensemble mean)
avg_model_year <- avg_model_year_long[, Mean:=mean(annual_sleep_loss_per_person), by=.(year)]
#avg_model_year

#plot annual projected sleep loss per person for 2050 and 2099
g<-ggplot(data=avg_model_year_long, aes(x=year, y=annual_sleep_loss_per_person, group=model)) +
  geom_line(color="purple", size=.5, alpha=0.3) + theme_classic() +
  scale_y_continuous(breaks=seq(40,65,5),limits=c(42.5,65.5))

v <- g + geom_line(data=avg_model_year, aes(x=year, y=Mean), color="purple", size= 1, alpha=1)

v
```

#Excess Sleep Loss RCP8.5
```{r}
#1.Compute and plot NASA SEDAC GPWv4 population-weighted avg. annual temperature-attributed EXCESS SLEEP LOSS for every model-by-year (2010,2050,2099) using NASA NEX BCSD CMIP5 (RCP8.5) spline sleep projections

#load sleep loss projection data for RCP8.5 (spline model prediction run on climate model time-series)
geocoords_labels<-fread(paste0(pathdir,"model_data/global_downscaled_sleeploss_projection_RCP8.5_popweights.csv"))
geocoords_labels[,V1:=NULL] #drop index col

#compute population weights using NASA SEDAC 2015 WP data
geocoords_labels[,pop_tot:=sum(pop_count)]
geocoords_labels[,pop_wts:=pop_count/pop_tot]
#geocoords_labels[,.(pop_count,pop_tot,pop_wts)]

#for each model, multiply sleep loss by pop weight fraction (note: need to specify climate model columns)
cols<-names(geocoords_labels)[3:65] #select model cols
#geocoords_labels
geocoords_labels[,(cols):=lapply(.SD, "*",pop_wts),.SDcols=cols]
#geocoords_labels

#sum weighted sleep loss values
avg_model_year <- geocoords_labels[, lapply(.SD, sum), .SDcols = cols]
#avg_model_year

#compute avg. for 2050 and 2099 across all vars
avg_model_year_long <- melt(avg_model_year)

#scale effect to sleep loss (in hours) per person
avg_model_year_long$transformed <- avg_model_year_long[,.(transformed = (value/60))]

#create year variable 
avg_model_year_long$year <- str_sub(avg_model_year_long$variable, -7, -4)

#create model variable 
avg_model_year_long$model <- str_sub(avg_model_year_long$variable, -17, -8)

#add 2010 baseline
avg_model_2010<-unique(avg_model_year_long[year=="2010",.(model,loss2010=transformed)])
avg_model_year_long[,.N]
avg_model_year_long<-merge(avg_model_year_long,avg_model_2010,all.x=T,by=c("model"))
avg_model_year_long[,.N]

setnames(avg_model_year_long, "transformed", "annual_sleep_loss_per_person") #rename y var

#compute total annual sleep loss by model-year
avg_model_year_long[year!=2010,annual_sleep_loss_per_person:=(loss2010 + annual_sleep_loss_per_person)]
#avg_model_year_long
```


```{r}
#inspect min, median and max sleep loss projection values by year
avg_model_year_long[year=="2010",.(med=median(annual_sleep_loss_per_person), low = min(annual_sleep_loss_per_person), high = max(annual_sleep_loss_per_person))]
avg_model_year_long[year=="2050",.(med=median(annual_sleep_loss_per_person), low = min(annual_sleep_loss_per_person), high = max(annual_sleep_loss_per_person))]
avg_model_year_long[year=="2099",.(med=median(annual_sleep_loss_per_person), low = min(annual_sleep_loss_per_person), high = max(annual_sleep_loss_per_person))]
```


```{r}
#RCP4.5 + RCP8.5 plot layers
#Plot population-weighted average excess sleep loss due to temperature (hours of sleep lost per person per year)
#average loss by year across all models (ensemble mean)
avg_model_year <- avg_model_year_long[, Mean:=mean(annual_sleep_loss_per_person), by=.(year)]
#avg_model_year

#plot annual projected excess sleep loss per person for 2050 and 2099
e<-ggplot(data=avg_model_year_long, aes(x=year, y=annual_sleep_loss_per_person, group=model)) +
  geom_line(color="orange", size=.5, alpha=0.3) + theme_classic() +
  scale_y_continuous(breaks=seq(40,65,5),limits=c(42.5,65.5))

f <- e + geom_line(data=avg_model_year, aes(x=year, y=Mean), color="orange", size= 1, alpha=1)

v
f 

#ggsave(paste0(pathdir,"model_data/annual_loss_pop_weighted_8.5.pdf"),width=4, height=4, dpi=300)
```

#Excess Short Sleep RCP4.5
```{r}
#1.Compute and plot NASA SEDAC GPWv4 population-weighted avg. annual temperature-attributed EXCESS Short <7hr. Sleep due to suboptimal temperature for every model-by-year (2010,2050,2099) using NASA NEX BCSD CMIP5 (RCP4.5) spline sleep projections

#load short sleep projection data for RCP4.5 (spline model prediction run on downscaled climate model time-series)
geocoords_labels<-fread(paste0(pathdir,"model_data/global_downscaled_shortsleep_projection_RCP4.5_popweights.csv"))
geocoords_labels[,V1:=NULL] #drop index col

#compute population weights using NASA SEDAC 2015 WP data
geocoords_labels[,pop_tot:=sum(pop_count)]
geocoords_labels[,pop_wts:=pop_count/pop_tot]
#geocoords_labels[,.(pop_count,pop_tot,pop_wts)]

#for each model, multiply sleep loss by pop weight fraction (note: need to specify climate model columns)
cols<-names(geocoords_labels)[3:65] #select model cols
#geocoords_labels
geocoords_labels[,(cols):=lapply(.SD, "*",pop_wts),.SDcols=cols]
#geocoords_labels

#sum weighted sleep loss values
avg_model_year <- geocoords_labels[, lapply(.SD, sum), .SDcols = cols]
#avg_model_year

#compute avg. for 2050 and 2099 across all vars
avg_model_year_long <- melt(avg_model_year)

#scale effect to absolute # of nights of short sleep per person
avg_model_year_long$transformed <- avg_model_year_long[,.(transformed = (value))]

#create year variable 
avg_model_year_long$year <- str_sub(avg_model_year_long$variable, -7, -4)

#create model variable 
avg_model_year_long$model <- str_sub(avg_model_year_long$variable, -17, -8)

#add 2010 baseline
avg_model_2010<-unique(avg_model_year_long[year=="2010",.(model,loss2010=transformed)])
avg_model_year_long[,.N]
avg_model_year_long<-merge(avg_model_year_long,avg_model_2010,all.x=T,by=c("model"))
avg_model_year_long[,.N]

setnames(avg_model_year_long, "transformed", "annual_short_sleep_per_person") #rename y var

#compute total annual # of nights of short sleep per person by model-year
avg_model_year_long[year!=2010,annual_short_sleep_per_person:=(loss2010 + annual_short_sleep_per_person)]
#avg_model_year_long
```



```{r}
#inspect min, median and max short sleep projection values by year
avg_model_year_long[year=="2010",.(med=median(annual_short_sleep_per_person), low = min(annual_short_sleep_per_person), high = max(annual_short_sleep_per_person))]
avg_model_year_long[year=="2050",.(med=median(annual_short_sleep_per_person), low = min(annual_short_sleep_per_person), high = max(annual_short_sleep_per_person))]
avg_model_year_long[year=="2099",.(med=median(annual_short_sleep_per_person), low = min(annual_short_sleep_per_person), high = max(annual_short_sleep_per_person))]
```

```{r}
#RCP4.5 layer
#Plot population-weighted average excess short sleep due to temperature (nights of short sleep per person per year)
#average short sleep by year across all models (ensemble mean)
 
gp<-ggplot(avg_model_year_long, aes(x = annual_short_sleep_per_person, y = year)) +
  geom_density_ridges(alpha=0.6, fill="purple") +
  theme_ridges() + 
  theme(legend.position = "none",axis.text = element_text(size = 9)) +
  scale_x_continuous(breaks=seq(10,17.5,.5),limits=c(10,18))

gp
```

#Excess Short Sleep RCP8.5
```{r}
#1.Compute and plot NASA SEDAC GPWv4 population-weighted avg. annual temperature-attributed EXCESS Short <7hr. Sleep due to suboptimal temperature for every model-by-year (2010,2050,2099) using NASA NEX BCSD CMIP5 (RCP8.5) spline sleep projections

#load short sleep projection data for RCP8.5 (spline model prediction run on climate model time-series)
geocoords_labels<-fread(paste0(pathdir,"/model_data/global_downscaled_shortsleep_projection_RCP8.5_popweights.csv"))
geocoords_labels[,V1:=NULL] #drop index col

#compute population weights using NASA SEDAC 2015 WP data
geocoords_labels[,pop_tot:=sum(pop_count)]
geocoords_labels[,pop_wts:=pop_count/pop_tot]
#geocoords_labels[,.(pop_count,pop_tot,pop_wts)]

#for each model, multiply sleep loss by pop weight fraction (note: need to specify climate model columns)
cols<-names(geocoords_labels)[3:65] #select model cols
#geocoords_labels
geocoords_labels[,(cols):=lapply(.SD, "*",pop_wts),.SDcols=cols]
#geocoords_labels

#sum weighted sleep loss values
avg_model_year <- geocoords_labels[, lapply(.SD, sum), .SDcols = cols]
#avg_model_year

#compute avg. for 2050 and 2099 across all vars
avg_model_year_long_2 <- melt(avg_model_year)

#scale effect to absolute # of nights of short sleep per person
avg_model_year_long_2$transformed <- avg_model_year_long_2[,.(transformed = (value))]

#create year variable 
avg_model_year_long_2$year <- str_sub(avg_model_year_long_2$variable, -7, -4)

#create model variable 
avg_model_year_long_2$model <- str_sub(avg_model_year_long_2$variable, -17, -8)

#add 2010 baseline
avg_model_2010<-unique(avg_model_year_long_2[year=="2010",.(model,loss2010=transformed)])
avg_model_year_long_2[,.N]
avg_model_year_long_2<-merge(avg_model_year_long_2,avg_model_2010,all.x=T,by=c("model"))
avg_model_year_long_2[,.N]

setnames(avg_model_year_long_2, "transformed", "annual_short_sleep_per_person") #rename y var

#compute total annual # of nights of short sleep per person by model-year
avg_model_year_long_2[year!=2010,annual_short_sleep_per_person:=(loss2010 + annual_short_sleep_per_person)]
#avg_model_year_long_2
```

```{r}
#inspect min, median and max short sleep projection values by year
avg_model_year_long_2[year=="2010",.(med=median(annual_short_sleep_per_person), low = min(annual_short_sleep_per_person), high = max(annual_short_sleep_per_person))]
avg_model_year_long_2[year=="2050",.(med=median(annual_short_sleep_per_person), low = min(annual_short_sleep_per_person), high = max(annual_short_sleep_per_person))]
avg_model_year_long_2[year=="2099",.(med=median(annual_short_sleep_per_person), low = min(annual_short_sleep_per_person), high = max(annual_short_sleep_per_person))]
```

```{r}
#RCP4.5 + RCP8.5 plot layers
#Plot population-weighted average excess short sleep due to temperature (nights of short sleep per person per year)
#average short sleep by year across all models (ensemble mean)

#combine RCP4.5 and 8.5 projections
projections<-list(avg_model_year_long,avg_model_year_long_2)
avg_model_year_long[,.N]
avg_model_year_long_2[,.N]
avg_model_year_long<-rbindlist(projections, use.names=T)
avg_model_year_long[,.N]

#create year-by-RCP variable
avg_model_year_long[,rcp:=substr(variable,17,21)]
avg_model_year_long[,year_rcp:=paste0(year,sep="_",substr(variable,17,21))]
#unique(avg_model_year_long$year_rcp)

gp<-ggplot(avg_model_year_long, aes(x = annual_short_sleep_per_person, y = year_rcp, fill = rcp, alpha=.5)) +
  geom_density_ridges(jittered_points = TRUE,
    position = position_points_jitter(width = 0.00, height = 0),
    point_shape = '|', point_size = 3, point_alpha = 1, alpha = 0.6) +
  stat_density_ridges(quantile_lines = TRUE, quantiles = 2) +
  theme_ridges() + 
  theme(legend.position = "none",axis.text = element_text(size = 9)) +
  scale_x_continuous(breaks=seq(10,17.5,.5),limits=c(10,18))

gp

#ggsave(paste0(pathdir,"model_data/annual_shortsleep_pop_weighted_8.5.pdf"),width=6, height=4, dpi=300)
```

##WORLD MAPS
#Net Sleep Loss RCP4.5 (net of 2010 loss)
```{r}
#World Map Sleep Loss RCP4.5
#1.Map avg. annual temperature-attributed NASA NEX BCSD CMIP5 (RCP4.5) ensemble mean NET SLEEP LOSS for each gridcell-by-year (2050-2010,2099-2010) using  spline sleep projections (for all countries in sleep data)

#load projection map data
sumdiff2050<-fread(paste0(pathdir,"model_data/net_sleep_loss_map_data_2050.csv"),header=T)
sumdiff2099<-fread(paste0(pathdir,"model_data/net_sleep_loss_map_data_2099.csv"),header=T)

#plot sleep projection maps
#set theme
my_theme <- theme_void()

#gradient scale color values
sleepalette2 <-c('#1B1349', '#CC33BA', '#EE6660','#FCC34F','#F9DF8C','#F9DF8C','#201B21')

#plot raster map
p1 <- ggplot() +
  geom_raster(data = sumdiff2050, aes(x = lon, y = lat, fill = Mean_2050)) +
scale_fill_gradientn(na.value = "grey50", values = scales::rescale(c(0, .2, .4, .6, .8, .999, .99999999, 1)), colours = sleepalette2, limits = c(0,26), position = 'left', breaks = c(seq(0, 26,  by=2))) + coord_sf(crs = st_crs(4326), xlim = c(-180, 180), ylim = c(-55, 80)) + my_theme

#add scale bar
p1 + guides(fill = guide_colourbar(barwidth = 32.9, barheight = 1, title = "Projected Annual Individual Sleep Loss Compared to 2010 Baseline (Hours Per Person)", title.position = "bottom", hjust = .5)) + theme(legend.position="bottom")

#ggsave(paste0(pathdir,"model_data/2050_World_Sleep_Loss_Map_RCP4.5.pdf"))

p2 <- ggplot() + geom_raster(data = sumdiff2099, aes(x = lon, y = lat, fill = Mean_2099)) + 
scale_fill_gradientn(na.value = "grey50", values = scales::rescale(c(0, .2, .4, .6, .8, .999,       .99999999, 1)), colours = sleepalette2, limits = c(0,26), position = 'left', breaks = c(seq(0, 26,  by=2))) + coord_sf(crs = st_crs(4326), xlim = c(-180, 180), ylim = c(-55, 80)) + my_theme

#add scale bar
p2 + guides(fill = guide_colourbar(barwidth = 32.9, barheight = 1, title = "Projected Annual Individual Sleep Loss Compared to 2010 Baseline (Hours Per Person)", title.position = "bottom", hjust = .5)) + theme(legend.position="bottom")

#ggsave(paste0(pathdir,"model_data/2099_World_Sleep_Loss_Map_RCP4.5.pdf"))
```


#Net Sleep Loss RCP8.5 (net of 2010 loss)
```{r}
#World Map Sleep Loss RCP8.5
#1.Map avg. annual temperature-attributed NASA NEX BCSD CMIP5 (RCP8.5) ensemble mean NET SLEEP LOSS for each gridcell-by-year (2050-2010,2099-2010) using  spline sleep projections (for all countries in data)

#load projection map data
sumdiff2050<-fread(paste0(pathdir,"model_data/net_sleep_loss_map_data_2050_RCP8.5.csv"),header=T)
sumdiff2099<-fread(paste0(pathdir,"model_data/net_sleep_loss_map_data_2099_RCP8.5.csv"),header=T)

#create sleep projection maps of pixel-level net sleep change due to climatic warming
#set theme
my_theme <- theme_void()

#gradient scale color values
sleepalette2 <-c('#1B1349', '#CC33BA', '#EE6660','#FCC34F','#F9DF8C','#F9DF8C','#201B21')

#plot raster map
p1 <- ggplot() +
  geom_raster(data = sumdiff2050, aes(x = lon, y = lat, fill = Mean_2050)) +
scale_fill_gradientn(na.value = "grey50", values = scales::rescale(c(0, .2, .4, .6, .8, .999, .99999999, 1)), colours = sleepalette2, limits = c(0,26), position = 'left', breaks = c(seq(0, 26,  by=2))) + coord_sf(crs = st_crs(4326), xlim = c(-180, 180), ylim = c(-55, 80)) + my_theme

#add scale bar
p1 + guides(fill = guide_colourbar(barwidth = 32.9, barheight = 1, title = "Projected Net Annual Individual Temperature-Attributed Sleep Loss - Net of 2010 Baseline (Hours Per Person Per Year) RCP8.5", title.position = "bottom", hjust = .5)) + theme(legend.position="bottom")

#ggsave(paste0(pathdir,"model_data/2050_World_Sleep_Loss_Map_RCP8.5.pdf"))

p2 <- ggplot() + geom_raster(data = sumdiff2099, aes(x = lon, y = lat, fill = Mean_2099)) + 
scale_fill_gradientn(na.value = "grey50", values = scales::rescale(c(0, .2, .4, .6, .8, .999, .99999999, 1)), colours = sleepalette2, limits = c(0,26), position = 'left', breaks = c(seq(0, 26,  by=2))) + coord_sf(crs = st_crs(4326), xlim = c(-180, 180), ylim = c(-55, 80)) + my_theme

#add scale bar
p2 + guides(fill = guide_colourbar(barwidth = 32.9, barheight = 1, title = "Projected Net Annual Individual Temperature-Attributed Sleep Loss - Net of 2010 Baseline (Hours Per Person Per Year) RCP8.5", title.position = "bottom", hjust = .5)) + theme(legend.position="bottom")

#ggsave(paste0(pathdir,"model_data/2099_World_Sleep_Loss_Map_RCP8.5.pdf"))
```


#Net Short Sleep RCP4.5 (net of 2010 short sleep)
```{r}
#World Map Short Sleep RCP4.5
#1.Map avg. annual temperature-attributed NASA NEX BCSD CMIP5 (RCP4.5) ensemble mean NET SHORT SLEEP for each gridcell-by-year (2050-2010,2099-2010) using  spline sleep projections (for all countries in sleep data)

#load projection map data
sumdiff2050<-fread(paste0(pathdir,"model_data/net_short_sleep_map_data_2050.csv"),header=T)
sumdiff2099<-fread(paste0("model_data/net_short_sleep_map_data_2099.csv"),header=T)

#create sleep projection maps of pixel-level net sleep change due to climate warming
#set theme
my_theme <- theme_void()

#gradient scale color values
sleepalette3 <-c('#0a0908',"#685e4d",'#ffe6b7','#f39a4f','#e76254','#b50909',"#840b14")

#plot raster map
p1 <- ggplot() +
  geom_raster(data = sumdiff2050, aes(x = lon, y = lat, fill = abs(Mean_2050))) +
scale_fill_gradientn(na.value = "grey50", values = scales::rescale(c(.000000001,.001,.2, .4, .6, .8, 1)), colours = sleepalette3, limits = c(0,7500), position = 'left', breaks = c(seq(0, 7000,  by=500))) + coord_sf(crs = st_crs(4326), xlim = c(-180, 180), ylim = c(-55, 80)) + my_theme

#add scale bar
p1 + guides(fill = guide_colourbar(barwidth = 32.9, barheight = 1, title = "Projected Net Annual Individual Temperature-Attributed Nights of Short <7hr. Sleep - Net of 2010 Baseline (Hours Per Person Per Year) RCP4.5", title.position = "bottom", hjust = .5)) + theme(legend.position="bottom")

##ggsave(paste0(pathdir,"model_data/2050_World_Short_Sleep_4.5_Map.pdf"))

p2 <- ggplot() + geom_raster(data = sumdiff2099, aes(x = lon, y = lat, fill = abs(Mean_2099))) +
scale_fill_gradientn(na.value = "grey50", values = scales::rescale(c(.000000001,.001,.2, .4, .6, .8, 1)), colours = sleepalette3, limits = c(0,7500), position = 'left', breaks = c(seq(0, 7000,  by=500))) + coord_sf(crs = st_crs(4326), xlim = c(-180, 180), ylim = c(-55, 80)) + my_theme

#add scale bar
p2 + guides(fill = guide_colourbar(barwidth = 32.9, barheight = 1, title = "Projected Net Annual Individual Temperature-Attributed Sleep Loss - Net of 2010 Baseline (Hours Per Person Per Year) RCP4.5", title.position = "bottom", hjust = .5)) + theme(legend.position="bottom")

#ggsave(paste0(pathdir,"model_data/2099_World_Short_Sleep_4.5_Map.pdf"))
```

#Net Short Sleep RCP8.5 (net of 2010 short sleep)
```{r}
#World Map Short Sleep RCP8.5
#1.Map avg. annual temperature-attributed NASA NEX BCSD CMIP5 (RCP4.5) ensemble mean NET SHORT SLEEP for each gridcell-by-year (2050-2010,2099-2010) using  spline sleep projections (for all countries in sleep data)

#load projection map data
sumdiff2050<-fread(paste0(pathdir,"model_data/net_short_sleep_map_data_2050_RCP8.5.csv"),header=T)
sumdiff2099<-fread(paste0(pathdir,"model_data/net_short_sleep_map_data_2099_RCP8.5.csv"),header=T)

#create sleep projection maps of pixel-level net sleep change due to climate warming
#set theme
my_theme <- theme_void()

#gradient scale color values
sleepalette3 <-c('#0a0908',"#685e4d",'#ffe6b7','#f39a4f','#e76254','#b50909','#840b14')

#plot raster map
p1 <- ggplot() +
  geom_raster(data = sumdiff2050, aes(x = lon, y = lat, fill = abs(Mean_2050))) +
scale_fill_gradientn(na.value = "grey50", values = scales::rescale(c(.000000001,.001,.2, .4, .6, .8, 1)), colours = sleepalette3, limits = c(0,7500), position = 'left', breaks = c(seq(0, 7000,  by=500))) + coord_sf(crs = st_crs(4326), xlim = c(-180, 180), ylim = c(-55, 80)) + my_theme

#add scale bar
p1 + guides(fill = guide_colourbar(barwidth = 32.9, barheight = 1, title = "Projected Net Annual Individual Temperature-Attributed Nights of Short <7hr. Sleep - Net of 2010 Baseline (Hours Per Person Per Year) RCP8.5", title.position = "bottom", hjust = .5)) + theme(legend.position="bottom")

#ggsave(paste0(pathdir,"model_data/2050_World_Short_Sleep_8.5_Map.pdf"))

p2 <- ggplot() + geom_raster(data = sumdiff2099, aes(x = lon, y = lat, fill = abs(Mean_2099))) +
scale_fill_gradientn(na.value = "grey50", values = scales::rescale(c(.000000001,.001,.2, .4, .6, .8, 1)), colours = sleepalette3, limits = c(0,7500), position = 'left', breaks = c(seq(0, 7000,  by=500))) + coord_sf(crs = st_crs(4326), xlim = c(-180, 180), ylim = c(-55, 80)) + my_theme

#add scale bar
p2 + guides(fill = guide_colourbar(barwidth = 32.9, barheight = 1, title = "Projected Net Annual Individual Temperature-Attributed Nights of Short <7hr. Sleep - Net of 2010 Baseline (Hours Per Person Per Year) RCP8.5", title.position = "bottom", hjust = .5)) + theme(legend.position="bottom")

#ggsave(paste0(pathdir,"model_data/2099_World_Short_Sleep_8.5_Map.pdf"))
```



